A model for the distribution of aftershock waiting times. 
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In this work the distribution of inter-occurrence times between earthquakes in aftershock sequences 
is analyzed and a model based on a non-homogeneous Poisson (NHP) process is proposed to quantify 
the observed scaling. In this model the generalized Omori's law for the decay of aftershocks is used as 
a time- dependent rate in the NHP process. The analytically derived distribution of inter-occurrence 
times is applied to several major aftershock sequences in California to confirm the validity of the 
proposed hypothesis. 



The occurrence of earthquakes is an outcome of com- 
plex nonlinear threshold dynamics in the brittle part of 
the Earth's crust. This dynamics is a combined effect of 
different temporal and spatial processes taking place in 
a highly heterogeneous media over a wide range of tem- 
poral and spatial scales 0. Despite this complexity, one 
can consider earthquakes as a point process in space and 
time, by neglecting the spatial scale of earthquake rup- 
ture zones and the temporal scale of the duration of each 
earthquake |^ |^ Then one can study the statisti- 
cal properties of this process and test models that may 
explain the observed seismic activity. 

In this letter, we analyze one of the important as- 
pects of the seismic activity, i.e. the statistics of inter- 
occurrences or waiting times of successive earthquakes 
in an aftershock sequence. To summarize our results, 
we have found that these statistics are consistent with a 
NHP process driven by a power-law decay rate by map- 
ping the occurrence of aftershocks in a multi-dimensional 
space into a marked point process on the one-dimensional 
time-line. The nature of this distribution is closely re- 
lated to the temporal correlations between earthquakes 
and can be used in hazard assessments of the occurrence 
of aftershocks. We have also derived an exact formula de- 
scribing the distribution of waiting times between events 
in a NHP process over a finite time period T and con- 
firmed it by numerical simulations. 

Earthquakes form a hierarchical structure in space and 
time and can also be thought of as a branching process 
where each event can trigger a sequence of secondary 
events and so forth. According to this structure, in some 
cases it is possible to discriminate between foreshocks, 
main shocks, and aftershocks, although, this classifica- 
tion is not well defined and can be ambiguous. However, 
it is observed that moderate and strong earthquakes ini- 
tiate sequences of secondary events which decay in time. 
These sequences are called aftershocks and their spatial 
and temporal distributions provide valuable information 
about the earthquake generating process (Sj]. 

Earthquakes follow several empirical scaling laws. 
Most prominent of them is Gutenberg-Richter relation 
which states that the cumulative number of events 



greater than magnitude m, A''(> m), follows an expo- 
nential distribution N{> m) (x 10"'"", where 6 is a uni- 
versal exponent near unity. This distribution becomes a 
power-law when the magnitude is replaced with the seis- 
mic moment, as the magnitude scales as a logarithm of 
the seismic moment. Another empirical power-law scal- 
ing relation describes the temporal decay of aftershock se- 
quences in terms of the frequency of occurrence of events 
per unit time, this is called the modified Omori's law 
1^. The spatial distribution of faults on which earth- 
quakes occur also satisfy (multi-)fractal statistics jSti^- 
These laws are manifestations of the self-similar nature 
of earthquake processes. 

Based on studies of properties of California seismic- 
ity, an attempt to introduce a unified scaling law for the 
temporal distribution of earthquakes was proposed |lCll| . 
The distribution of inter-occurrence times between suc- 
cessive earthquakes was obtained by using as scaling pa- 
rameter both a grid size over which the region was subdi- 
vided, and a lower magnitude cutoff. Two distinct scaling 
regimes were found. For short times, aftershocks domi- 
nate the scaling properties of the distribution, decaying 
according to the modified Omori's law. For long times, 
an exponential scaling regime was found that can be as- 
sociated with the occurrence of main shocks. To take 
into account the spatial heterogeneity of seismic activity, 
it has been argued that the second regime is not an ex- 
ponential but another power-law [jjj . An analysis of the 
change in behavior between these two regimes based on a 
nonstationary Poisson sequence of events was carried out 
in 12]. The further analysis of aftershock sequences in 
California and Iceland revealed the existence of another 
scaling regime for small values of inter-occurrence times 



An alternative approach to describe a unified scaling of 
earthquake recurrence times was suggested in [l^[l3 . ll5j . 
where the distributions computed for different spatial ar- 
eas and magnitude ranges were rescaled with the rate 
of seismic activity for each area considered. It was ar- 
gued that the seismic rate fully controls the observed 
scaling, and that the shape of the distribution can be ap- 
proximated by the generalized gamma function, indicat- 
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ing the existence of correlations in the recurrence times 
beyond the duration of aftershock sequences. This ap- 
proach agrees with observations that main shocks show 
a nonrandom behavior with some effects of long-range 
memory [l^ . 

Before we carry out an analysis of seismic data, we will 
outline a derivation of a distribution function for waiting 
times between events in a point process characterized by 
a rate r{t) and distributed according to NHP statistics 
over a finite time interval T. The full analysis is going 
to be reported elsewhere 0. The instantaneous proba- 
bility distribution function of waiting times Y at time t, 
until the next event in accordance with the NHP process 
hypothesis has the following form |3| 



F{t, At) = Prob{y < At} = 1 - c /t 



r(u) du 



(1) 



where r(t) is a rate of occurrence of events at the time 
t. The probability density function of waiting times over 
the whole time period T has the following form 
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where N — r(u) du is the total number of events 
during a time period T. In the simple case of a 
constant rate (r = const) one recovers the result for 
the homogeneous Poisson process, namely, PriAt) = 



-r At 



rAt 
T 
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In order to check the correctness of the derivation of 
Eq. |(5J) we have performed numerical simulations of the 
NHP process with a decaying time dependent rate r{t). 
Specifically, we have used a power-law rate defined as 



rit) = 



1 



r (l + t/cf ' 



(3) 



where r is a characteristic time that defines the rate at 
time t — 0, c is a second characteristic time that elim- 
inates the singularity at t = 0, and p is a power-law 
exponent. This rate is commonly used to describe the 
relaxation of aftershock sequences after a main shock and 
is called the modified Omori's law Q. 

In Fig. ^ we show plots of numerical simulations of 
the NHP process with varying scaling parameters t, c, 
and T with fixed p = 1.2. These are indicated as solid 
symbols. We also plot the corresponding numerical in- 
tegrations of Eq. (0) for the same values of r, c, and T. 
The comparison shows that Eq. jSJ correctly describes 
the inter-occurrence time distributions between events. 

We have also performed a scaling analysis of our sim- 
ulated distributions for the values At > 150 s. This is 
shown in Fig. ^p. The distributions collapse onto each 
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c = 1 day, T = 45 days 
= 2 s, c = 1 day, T = 360 days 
= 10 s, c = 1 day, T = 180 days 
= 2 s, c = 1 /1 6 day, T = 90 days 
= 2 s, c= 1/4 day, T = 720 days 
= 50 s, c = 1/4 day, T = 90 days 
= 1/4day, T = 720 days 
= 1 day, T = 1440 days 
= 1/16day,T = 720days 



r > T = 10s,c = 1/16day,T=1440days 




FIG. 1: a) Inter-occurrence time distributions for the NHP 
process with the decaying rate (j^J. Symbols are numerical 
simulations of the NHP process. Solid lines are numerical 
integrations of Eq. Q. b) The scaling analysis of the distri- 
butions for At > 150 s using Eq. Q. The solid line is a fit of 
the generalized gamma function to the rescaled data. 



other with respect to the following scaling law 



1 /At 



^At) = - — 



At /T 
r V c 
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(4) 



where a = 1.212 and (3 — 1.194. The scaling function 
f{x) can be approximated by the generalized gamma 
function f{x) — Ax'' exp{~x/ B) with 7 — 0.056, A — 
0.329, and B = 1.122. The distribution functions have 
two characteristic time-scales, i.e. r and t{T/c)^, which 
define two roll overs for small and large At's. 

To find scaling relations between the above exponents, 
one has to perform an asymptotic analysis of the integral 
(0) for large At ^ 1 and T 1. This analysis can be 
done by using Laplace's method and gives P{At) ~ 
T ^ 00 and P{At) ~ exp[-At/T(l + 
T/c)'P] for finite T, from which we conclude that a — 7 = 
2—1/p and [3 = p. The exponent 2— 1/p was also reported 
in J. 

To check the proposed hypothesis that aftershocks can 
be modeled as a NHP process, the derived distribution 
^ has been compared to several aftershock sequences in 
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California. We have used the seismic catalog provided 
by the Southern California Earthquake Center (SCSN 
catalog, http://www.data.scec.org). The identification 
of aftershock sequences is usually an ambiguous problem 
which requires assumptions on the spatial and tempo- 
ral clustering of aftershocks Q . In this work we assume 
that all events that have occurred after a main shock 
within a given time interval T and a square area L x L 
are aftershocks of a particular main shock. By neglecting 
the magnitude, spatial size and duration of each individ- 
ual event and considering all aftershocks above a certain 
threshold rric, we map a multi-dimensional process into 
a process on the one-dimensional time-line. This marked 
process is characterized by the times of occurrence of in- 
dividual events U and their magnitudes rrii. We define 
inter-occurrence or waiting times between successive af- 
tershocks as Ati = ti — ti-i, 1 = 1,2,... and study their 
statistical distribution over a finite time interval T. 

In this work, we use the decay rate of aftershocks, 
introduced in Eoj , which is a generalization of the 
modified Omori's law 0, where the characteristic time 
c{m) — T {p — 1)10''^™ ""^^ is a function of the lower 
cutoff magnitude of an aftershock sequence m, m* is the 
maximum value of an aftershock in a sequence with finite 
number of events, and p > 1. We also assume that af- 
tershock sequences satisfy the Gutenberg-Richter cumu- 
lative distribution N{> m) = iqM™*-™). This defines 
a truncated distribution where the expected number of 
events with magnitudes greater than m* is eq ual to one 
and the exponent b is generally near unity |2lj . 

In Fig. 13 we have computed the distribution functions 
of inter-occurrence times between successive aftershocks 
from the observed data of three California aftershock se- 
quences. These are indicated as solid symbols in Fig. El 
For each of these sequences we have used square boxes 
of size L X L = 1.25° x 1.25° for the Landers earth- 
quake {m„,s = 7.3; Jun 28, 1992) and 1.0° x 1.0° for 
the Northridge {m„is — 6.7; Jan 17, 1994) and Hector 
Mine {m,ns = 7.1; Oct 16, 1999) earthquakes centered 
on the epicenter of the main shocks and a time interval 
of T = 1 year. All the earthquakes that occurred in the 
spatio-temporal boxes were treated as aftershocks. The 
analysis of the data shows that the distributions are not 
too sensitive to changes in the linear size L of the box, 
the results are almost the same for L ranging from 0.25° 
to 1.75° within statistical errors. This means that the 
distributions are dominated by the activity of the events 
generated by the main shocks and the background seis- 
micity doesn't contribute significantly to the scaling. 

In this analysis of an aftershock sequence as a point 
process we treat all earthquakes as having the same mag- 
nitude and as a result we lose a significant fraction of 
information related to the physics of the process. To re- 
cover some information from the magnitude domain of 
each sequence we have used a lower magnitude cutoff rric 
as a scaling parameter and study sequences with differ- 
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FIG. 2: Inter-occurrence time distributions for the three 
California aftershock sequences for different magnitude cut- 
offs nic = 1.5, 2.0, 2.5, 4.0. Square areas considered are 
1.25° X 1.25° for (a) the Landers sequence and 1.0° x 1.0° 
for (b) the Northridge and (c) Hector Mine sequences each 
centered on the epicenters of the main shocks. In all cases a 
time interval of T = 1 year following the main shock has been 
used. The solid lines have been computed using Eq. ||5J. 



ent Trie's. These are depicted by different colors in the 
plots (Fig. 12). The distributions with lower magnitude 
cutoffs have a shorter power-law scaling regime and start 
to roll over more quickly. This can be explained by the 
presence of more events in the sequences with lower mag- 
nitude cutoffs and as a result the shortening of the mean 
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□ Landers 
O Northridge 
A Hector Mine 



type of distribution that governs the occurrence of after- 
shocks is important in any hazard assessment programs. 
The derived distribution (Eq. has much broader ap- 
pUcabiUty and can be used for studies of many time de- 
pendent processes which follow NHP statistics. 
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FIG. 3: The scaling analysis of the inter-occurrence time dis- 
tributions of the three aftershocks sequences according to 
The red, green, and blue symbols correspond to the Landers 
(p = 1.22, T = 107.25 s), Northridge [p = 1.17, r = 53.14 s), 
and Hector Mine (p = 1.22, r — 83.71 s) sequences. In each 
sequence the varying parameters c and T = 180, 360, and 720 
days have been used. The best collapse have been found for 
a = 1.25 and f3 = 1.22. 

time intervals between events. Another scaling parame- 
ter which affects the roll over is the time interval T. An 
increase in T leads to the occurrence of longer intervals 
At between events. 

To compare the observed scaling with the simulations 
of a NHP process we also plot in Fig. |51 as solid curves, 
the distributions computed assuming that aftershock se- 
quences follow NHP statistics with the rate given by 
Eq. Q and the parameters r, c, and p estimated from the 
observed three California aftershock sequences [l^ . The 
plots show that the modeled distributions are in excellent 
agreement with the observations. 

We have also performed a scaling analysis of the 
inter-occurrence time distributions of these aftershock se- 
quences. This is shown in Fig|21 The characteristic time 
c and the time interval T have been chosen as scaling pa- 
rameters. These sequences are characterized by slightly 
different initial rates r and exponents p. The results show 
a reasonably good scaling with respect to c and T which 
supports our hypothesis that aftershock sequences can be 
described as a NHP process. 

In summary, the studies of inter-occurrence of after- 
shocks presented in this work suggest that aftershock 
sequences can be modeled to a good approximation as 
a point process governed by NHP statistics, where the 
rate of activity decays as a power-law (Eq. O . This de- 
caying rate introduces a self-similar regime into the ob- 
served scaling followed by an exponential roll over. An 
analysis of a nonstationary earthquake sequence was also 
performed in J^]. It was suggested the existence of a 
secondary clustering structure within the main sequence 
and deviation from NHP behavior. The knowledge of the 
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